!> author: 左志华
!> date: 2022-07-16
!>
!> 测试密度求和
module test_sph_density_summation

    use testdrive, only: new_unittest, unittest_type, error_type, check
    use sph_kinds, only: rk
    use sph_region, only: region_t
    use sph_density_summation, only: denormalized_density_summation, normalized_density_summation
    use twin_sph_list_type, only: twin_list_t
    use sph_list_type, only: list_t
    use sph_linked_list_nnps, only: build => build_grid_from_loc, query
    use sph_smoothed_kernel_function, only: skf
    use sph_env, only: nums
    implicit none
    private

    public :: collect_density_summation

contains

    subroutine collect_density_summation(testsuite)
        type(unittest_type), allocatable, intent(out) :: testsuite(:)

        allocate (testsuite, source=[ &
                  new_unittest("FUNC: denormalized_density_summation", test_density_summation), &
                  new_unittest("FUNC: normalized_density_summation", test_regularize_density_summation) &
                  ])

    end subroutine collect_density_summation

    subroutine test_density_summation(error)
        type(error_type), intent(out), allocatable :: error
        type(region_t) :: region
        type(list_t), allocatable :: grids(:, :)
        type(twin_list_t) :: pairs
        real(rk) :: hv(2), rho_expected(9)
        integer :: i

        ! 构建粒子域, 背景网格间距为 1，光滑长度为 1，粒子数为 9
        ! 3 6 9
        ! 2 5 8
        ! 1 4 7
        nums(1) = 9
        nums(4) = 9
        allocate (region%particles%loc(2, 9), region%particles%mass(9), region%particles%rho(9), region%particles%itype(9))
        region%particles%hsml = 1.0_rk
        region%particles%loc(1:2, 1) = [0.5_rk, 0.5_rk]
        region%particles%loc(1:2, 2) = [0.5_rk, 1.5_rk]
        region%particles%loc(1:2, 3) = [0.5_rk, 2.5_rk]
        region%particles%loc(1:2, 4) = [1.5_rk, 0.5_rk]
        region%particles%loc(1:2, 5) = [1.5_rk, 1.5_rk]
        region%particles%loc(1:2, 6) = [1.5_rk, 2.5_rk]
        region%particles%loc(1:2, 7) = [2.5_rk, 0.5_rk]
        region%particles%loc(1:2, 8) = [2.5_rk, 1.5_rk]
        region%particles%loc(1:2, 9) = [2.5_rk, 2.5_rk]
        region%particles%mass = 1000.0_rk
        region%particles%itype = 2

        ! 构建粒子对
        call build(region%particles%loc, 1.0_rk, 2, grids) ! scale 取为 1，因为 L 型查找最低网格要求是 3*3
        call query(region%particles%loc, grids, 2, 100, pairs)
        skf%type = 'cubic_spline'
        call skf%setup_kernel(region%particles%hsml)
        call denormalized_density_summation(region%particles, pairs, nums(1))

        call check(error, pairs%log, 20, more="pairs%log")  ! 粒子对为 20 对
        if (allocated(error)) return

        rho_expected = [704.9_rk, 841.5_rk, 704.9_rk, 841.5_rk, 1118.7854278159241_rk, 841.5_rk, 704.9_rk, &
                        841.5_rk, 704.9_rk]
        do i = 1, 9
            call check(error, region%particles%rho(i), rho_expected(i), &
                       more="region%particles%rho", thr=0.1_rk, rel=.true.)
            if (allocated(error)) return
        end do

    end subroutine test_density_summation

    subroutine test_regularize_density_summation(error)
        type(error_type), intent(out), allocatable :: error
        type(region_t) :: region
        type(list_t), allocatable :: grids(:, :)
        type(twin_list_t) :: pairs
        real(rk) :: hv(2), rho_expected(9)
        integer :: i

        ! 构建粒子域, 背景网格间距为 1，光滑长度为 1，粒子数为 9
        nums(1) = 9
        nums(4) = 9
        allocate (region%particles%loc(2, 9), region%particles%mass(9), region%particles%rho(9), region%particles%itype(9))
        region%particles%hsml = 1.0_rk
        region%particles%loc(1:2, 1) = [0.5_rk, 0.5_rk]
        region%particles%loc(1:2, 2) = [0.5_rk, 1.5_rk]
        region%particles%loc(1:2, 3) = [0.5_rk, 2.5_rk]
        region%particles%loc(1:2, 4) = [1.5_rk, 0.5_rk]
        region%particles%loc(1:2, 5) = [1.5_rk, 1.5_rk]
        region%particles%loc(1:2, 6) = [1.5_rk, 2.5_rk]
        region%particles%loc(1:2, 7) = [2.5_rk, 0.5_rk]
        region%particles%loc(1:2, 8) = [2.5_rk, 1.5_rk]
        region%particles%loc(1:2, 9) = [2.5_rk, 2.5_rk]
        region%particles%mass = 1000.0_rk
        region%particles%rho = 1000.0_rk
        region%particles%itype = 2

        ! 构建粒子对
        call build(region%particles%loc, 1.0_rk, 2, grids)
        call query(region%particles%loc, grids, 2, 100, pairs)
        skf%type = 'cubic_spline'
        call skf%setup_kernel(region%particles%hsml)
        call normalized_density_summation(region%particles, pairs, nums(1))

        call check(error, pairs%log, 20, more="pairs%log")
        if (allocated(error)) return

        rho_expected = [1000.00006_rk, 1000.00012_rk, 1000.00012_rk, 1000.00006_rk, &
                    &1000.00012_rk, 1000.00012_rk, 1000.00006_rk, 1000.00018_rk, 1000.00012_rk]
        do i = 1, 9
            call check(error, region%particles%rho(i), rho_expected(i), &
                       more="region%particles%rho", thr=0.1_rk, rel=.true.)
            if (allocated(error)) return
        end do

    end subroutine test_regularize_density_summation

end module test_sph_density_summation
